/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of LINAL (LINear ALgebra)
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "linal/common.hxx"
#include "linal/util.hxx"

#include "linal/format.hxx"

#include "linal/operation.hxx"

#include "linal/device/host/operation.hxx"

namespace linal {
  int qr(Hier_ A, Hier_ T) {
    if (!A.get_m() || !A.get_n()) return true;
    
    FLA_Obj ATL,   ATR,      A00, A01, A02, 
            ABL,   ABR,      A10, A11, A12,
                             A20, A21, A22;

    FLA_Obj TTL,   TTR,      T00, T01, T02, 
            TBL,   TBR,      T10, T11, W12,
                             T20, T21, T22;

    FLA_Part_2x2( ~A, &ATL, &ATR,
                      &ABL, &ABR, 0, 0, FLA_TL );

    FLA_Part_2x2( ~T, &TTL, &TTR,
                      &TBL, &TBR, 0, 0, FLA_TL );

    while ( FLA_Obj_min_dim( ABR ) > 0 ){
      FLA_Repart_2x2_to_3x3( ATL,  ATR,       &A00,  &A01, &A02,
                                              &A10,  &A11, &A12,
                             ABL,  ABR,       &A20,  &A21, &A22,
                             1, 1, FLA_BR );

      FLA_Repart_2x2_to_3x3( TTL,  TTR,       &T00,  &T01, &T02,
                                              &T10,  &T11, &W12,
                             TBL,  TBR,       &T20,  &T21, &T22,
                             1, 1, FLA_BR );

      //------------------------------------------------------------
      Hier_ 
        A_10(A10), A_11(A11), A_20(A20), A_21(A21), 
        T_01(T01), T_11(T11), W_12(W12); 

      FLA_Obj AB1;
      FLA_Merge_2x1( A11,
                     A21, &AB1 );

      Hier_ AB_1(AB1);
      FLA_QR_UT( ~AB_1.flat(), ~T_11.flat() );

      if ( FLA_Obj_width( A12 ) > 0 )  {
        FLA_Obj AB2;
        FLA_Merge_2x1( A12,
                       A22, &AB2 );

        Hier_ AB_2(AB2);

#pragma omp task
        {
          apply_q( FLA_LEFT, FLA_CONJ_TRANSPOSE, 
                   FLA_FORWARD, FLA_COLUMNWISE,
                   AB_1, T_11, W_12, AB_2);
        }

      }

      if ( A_10.get_m() && A_10.get_n() ) 
        FLA_Copyt( FLA_CONJ_TRANSPOSE, ~A_10.flat(), ~T_01.flat() );

#pragma omp task
      {
        trmm( FLA_RIGHT, FLA_LOWER_TRIANGULAR,
              FLA_NO_TRANSPOSE, FLA_UNIT_DIAG,
              FLA_ONE, A_11, T_01 );
          
        gemm( FLA_CONJ_TRANSPOSE, FLA_NO_TRANSPOSE, 
              FLA_ONE, A_20, A_21, FLA_ONE, T_01 );
      }
        
#pragma omp taskwait
	
      //------------------------------------------------------------
      FLA_Cont_with_3x3_to_2x2( &ATL,  &ATR,       A00, A01,  A02,
                                                   A10, A11,  A12,
                                &ABL,  &ABR,       A20, A21,  A22,
                                FLA_TL );
      FLA_Cont_with_3x3_to_2x2( &TTL,  &TTR,       T00, T01,  T02,
                                                   T10, T11,  W12,
                                &TBL,  &TBR,       T20, T21,  T22,
                                FLA_TL );
    }

    return true;
  }
}
